from sympy import *
from sympy.matrices import Matrix
from sympy.physics.vector import ReferenceFrame, curl, divergence, gradient
import plotly.graph_objects as go
import numpy as np
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt

Helper Functions

Curl and grad helper functions.

This code still needs testing functions written for it.

def nabla(tuplefield, kind='curl'):
    '''
    tuplefield: tuple
        tuple of expressions involving x,y,z
    '''
    R = ReferenceFrame('R')
    try:
        x,y,z = sorted(list(tuplefield.free_symbols), key=lambda i: i.name)
    except:
        return NameError('x,y,z not defined')
    # turn tuple into vector field using reference frame
    field = rhopf.subs({x:R[0],y:R[1],z:R[2]})
    field = field[0] * R.x + field[1] * R.y + field[2] * R.z
    if kind == 'curl':
        field = curl(field, R)
    if kind == 'divergence':
        field = divergence(field, R)
    # turn vector field back into vector
    field = field.subs({R[0]:x,R[1]:y,R[2]:z})
    field = (field.dot(R.x), field.dot(R.y), field.dot(R.z))
    return field

Stereographic Projection

This is a class that performs a stereographic projection between $\mathbb{R}^n$ and $S^n$.

class StereographicProjection:
    def __init__(self, special_point = (0,0,0,1), sphere_centre=(0,0,0,0), sphere_radius=1, n=3):
        '''
        special_point: tuple, default: (0,0,1)
            special reference point on Sn from which to perform stereographic projection
        sphere_centre: tuple, default: (0,0,0)
            centre of the sphere in R^{n+1}
        sphere_radius: float, default: 1
            radius of the sphere
        n: int, default: 2
            dimension of the plane of stereographic projection
        '''
        self.sp = Matrix(special_point)
        self.sc = Matrix(sphere_centre)
        self.r = sphere_radius
        self.n = n
        
        # assert special point lies on sphere
        assert (self.sp - self.sc).dot(self.sp - self.sc) == self.r
        
        # assert special point does not lie on Rn
        assert self.sp[-1] != 0
        
    def project_Rn_to_Sn(self, v):
        '''
        vector on Rn expressed as an (n+1)-tuple with last entry null
        '''
        v = Matrix(v)
        
        # assert point lies on Rn
        assert len(v) == self.n + 1
        assert v[-1] == 0
        
        t = Symbol('t')
        w = self.sp + t * (v - self.sp) - self.sc
        t = solve(w.dot(w) - self.r, t)
        if not t: return nan
        return tuple(simplify(self.sp + t[1] * (v - self.sp)))
    
    def project_Sn_to_Rn(self, v):
        '''
        vector on Sn expressed as an (n+1)-tuple
        '''
        v = Matrix(v)
        
        # assert point lies on Sn
        test = N(v.dot(v))
        if not test.free_symbols:
            assert np.isclose(float(test), 1)

        t = self.sp[-1] / (self.sp[-1] - v[-1])
        v = self.sp + t * (v - self.sp)
        
        if not test.free_symbols:
            return tuple([float(i) for i in N(v)])
        return tuple(simplify(v))

Riemann Hopf Map

This is the version of the Hopf map that is formulated using the Riemann sphere.

The Riemann_Hopf_Map function returns a vector field by performing the following series of maps:

$\mathbb{R}^3 \xrightarrow{\text{stereo}} S^3 \simeq \mathbb{C}^2 \xrightarrow{z_1/z_2} \mathbb{C} \simeq \mathbb{R}^2 \xrightarrow{\text{stereo}} S^2$.

def Riemann_Hopf_Map():
    # consider any point (x,y,z) in R3
    x,y,z = symbols('x y z', real=True)
    # stereographic projection to S3
    stereo1 = StereographicProjection(special_point = (0,0,0,1), sphere_centre=(0,0,0,0), sphere_radius=1, n=3)
    a,b,c,d = stereo1.project_Rn_to_Sn((x,y,z,0))
    # compute the Riemann Hopf map by dividing two complex numbers
    w = (a + I*b) / (c + I*d)
    rew = simplify(re(w))
    imw = simplify(im(w)) 
    # stereographic projection to S2
    stereo2 = StereographicProjection(special_point=(0,0,1), sphere_centre=(0,0,0), sphere_radius=1, n=2)
    return factor(stereo2.project_Rn_to_Sn((rew,imw,0)))

rhopf = Riemann_Hopf_Map(); rhopf

$\displaystyle \left( \frac{4 \left(x^{2} y + 2 x z + y^{3} + y z^{2} - y\right)}{\left(x^{2} + y^{2} + z^{2} + 1\right)^{2}}, \ - \frac{4 \left(x^{3} + x y^{2} + x z^{2} - x - 2 y z\right)}{\left(x^{2} + y^{2} + z^{2} + 1\right)^{2}}, \ - \frac{x^{4} + 2 x^{2} y^{2} + 2 x^{2} z^{2} - 6 x^{2} + y^{4} + 2 y^{2} z^{2} - 6 y^{2} + z^{4} + 2 z^{2} + 1}{\left(x^{2} + y^{2} + z^{2} + 1\right)^{2}}\right)$

Preimages

The preimage of a point under the Hopf map is a closed loop. The preimages corresponding to any two distinct points are a pair of linked closed loops.

We can explicitly compute the preimages of the Hopf map by noticing $$\eta(z_1,z_2)=\dfrac{z_1}{z_2} = \dfrac{r_1}{r_2}e^{i(\theta_1-\theta_2)} = R e^{i\phi}$$ And since $r_1^2 + r_2^2 = 1$, we have that $$z_1 = \sqrt{\dfrac{1}{R^2+1}}e^{i(\phi+t)},\quad z_2 = \sqrt{\dfrac{R^2}{R^2+1}}e^{it},$$ where $t$ is an arbitrary real parameter. Clearly $(z_1, z_2)$ is $2\pi$ periodic in $t$ and therefore forms a closed loop.

def Inverse_Riemann_Hopf_Map(vector):
    '''
    Given any tuple (x,y,z) in S2, computes the preimage in R3 under the Riemann Hopf Map.
    '''
    
    # map S2 to R2 equivalent to C
    stereo1 = StereographicProjection(special_point=(0,0,1), sphere_centre=(0,0,0), sphere_radius=1, n=2)
    a, b, _ = stereo1.project_Sn_to_Rn(vector)
    
    # compute preimages of Hopf map in C2
    
    # Rewrite a + bi in polar form R * exp(i * Phi)
    R = sqrt(a**2 + b**2)
    if N(a).free_symbols or N(b).free_symbols:
        Phi = atan(b/a)
    elif np.isclose(float(abs(b / a)), float(pi / 2)):
        Phi = [-1, 1][b / a >= 0] * pi / 2
    else:
        Phi = atan(b/a)
    # Compute radii of preimages
    r_1 = R / (R**2 + 1) ** .5
    r_2 = 1 / (R**2 + 1) ** .5
    # Write parametrised preimages
    t = Symbol('t')
    vector = (r_1 * cos(Phi + t), r_1 * sin(Phi + t), r_2 * cos(t), r_2 * sin(t))
    
    # map C2 to R3
    stereo2 = StereographicProjection(special_point=(0,0,0,1), sphere_centre=(0,0,0,0), sphere_radius=1, n=3)
    a, b, c, _ = stereo2.project_Sn_to_Rn(vector)
    
    return (a, b, c)

def compute_hopf_preimages(vector):
    '''
    Compute a set of Hopf map preimage points for a tuple in R3.
    '''
    vector = Inverse_Riemann_Hopf_Map(vector)

    t, theta = symbols('t theta')

    vectorx = lambdify((t, theta), vector[0], 'numpy')
    vectory = lambdify((t, theta), vector[1], 'numpy')
    vectorz = lambdify((t, theta), vector[2], 'numpy')

    mesht = np.mgrid[0:2*np.pi:300j]
    meshtheta = np.mgrid[0:2*np.pi:300j]

    x = vectorx(mesht, meshtheta)
    y = vectory(mesht, meshtheta)
    z = vectorz(mesht, meshtheta)
    
    return (x,y,z)

x, y, z = compute_hopf_preimages((1,0,0))
a, b, c = compute_hopf_preimages((sqrt(2)/2, 0, sqrt(2)/2))

fig = go.Figure()

fig.add_trace(go.Scatter3d(
    x=x, y=y, z=z,
    marker=dict(
        size=0,
    ),
    line=dict(
        color='darkblue',
        width=2
    )
))

fig.add_trace(go.Scatter3d(
    x=a, y=b, z=c,
    marker=dict(
        size=0,
    ),
    line=dict(
        color='darkred',
        width=2
    )
))

fig.update_layout(
    template='ggplot2',
    scene = dict(
        xaxis = dict(
            showbackground=False,
            showticklabels=False,
            showaxeslabels=False,
             ),
        yaxis = dict(
            showbackground=False,
            showticklabels=False,
            showaxeslabels=False,
            ),
        zaxis = dict(
            showbackground=False,
            showticklabels=False,
            ),
    ),
    width=700,
    margin=dict(
    r=10, l=10,
    b=10, t=10)
)

fig.show()

Hopf Invariant

The linking number of the preimages of any two distinct points for the standard Hopf map is a conserved quantity. This is formalised through the the Hopf invariant which can be defined as $$H(f) = \int_{S^3}\omega\wedge d\omega$$ where $\alpha$ is a generator of $H_{\text{DR}}^2(S^3)$ and $\omega$ is a $1$-form such that $f^*\alpha=d\omega$ and $f$ is a map from $S^3$ to $S^2$.

Programmatic implementation of this integral coming soon...

Hopf Texture

The Hopf texture can be visualised using the PT construction. We set a plane of orientation, say $n_z=0$, and plot the surface corresponding to vectors pointing along this orientation coloured by their phase.

def compute_hopf_texture():
    '''
    Draws the PT construction for the Riemann Hopf Map.
    '''
    # declare coordinate symbols and parametric symbols
    x,y,z = symbols('x y z')
    t, theta = symbols('t theta')

    # compute the preimages of the z = 0 orientation
    vector = Matrix(Inverse_Riemann_Hopf_Map((x,y,0))).subs({x**2+y**2:1, atan(y/x):theta})

    vectorx = lambdify((t, theta), vector[0], 'numpy')
    vectory = lambdify((t, theta), vector[1], 'numpy')
    vectorz = lambdify((t, theta), vector[2], 'numpy')
    
    mesht, meshtheta = np.mgrid[0:2*np.pi:100j, 0:2*np.pi:100j]

    x = vectorx(mesht, meshtheta)
    y = vectory(mesht, meshtheta)
    z = vectorz(mesht, meshtheta)
    
    # alongside coordinates return the orientation modulo pi as a colour scale
    return x, y, z, 2 * (meshtheta % np.pi)

x, y, z, colour = compute_hopf_texture()

fig = go.Figure(
    go.Surface(x=x, y=y, z=z, surfacecolor=colour, colorscale='greys'),
)

fig.update_layout(
    template='ggplot2',
    scene = dict(
        xaxis = dict(
            showbackground=False,
            showticklabels=False,
            showaxeslabels=False,
             ),
        yaxis = dict(
            showbackground=False,
            showticklabels=False,
            showaxeslabels=False,
            ),
        zaxis = dict(
            showbackground=False,
            showticklabels=False,
            ),
    ),
    width=700,
    margin=dict(
    r=10, l=10,
    b=10, t=10)
)

fig.show()

Twist Isosurface

The corresponding twist, $\mathbf{n} \cdot \nabla \times \mathbf{n}$, can be visualised through its isosurfaces.

twist = factor(Matrix(nabla(Matrix(rhopf), kind='curl')).dot(Matrix(rhopf))); twist

$\displaystyle \frac{8 \left(x^{2} z - 3 x^{2} + y^{2} z - 3 y^{2} + z^{3} + z^{2} + z + 1\right) \left(x^{2} z + 3 x^{2} + y^{2} z + 3 y^{2} + z^{3} - z^{2} + z - 1\right)}{\left(x^{2} + y^{2} + z^{2} + 1\right)^{4}}$

x,y,z = symbols('x y z')

numpy_twist = lambdify((x,y,z), twist, 'numpy')

X, Y, Z = np.mgrid[-5:5:40j, -5:5:40j, -5:5:40j]

values = numpy_twist(X,Y,Z)

fig = go.Figure(
    data=go.Isosurface(
    x=X.flatten(),
    y=Y.flatten(),
    z=Z.flatten(),
    value=values.flatten(),
    opacity=0.9,
    surface_fill=0.4,
    isomin=-1,
    isomax=1,
    surface_count=3,
    caps=dict(x_show=False, y_show=False),
    colorscale='Cividis',
    ),
)

fig.update_layout(
    scene = dict(
        xaxis = dict(
            showbackground=False,
            showticklabels=False,
            showaxeslabels=False,
             ),
        yaxis = dict(
            showbackground=False,
            showticklabels=False,
            showaxeslabels=False,
            ),
        zaxis = dict(
            showbackground=False,
            showticklabels=False,
            ),
    ),
    width=700,
    margin=dict(
    r=10, l=10,
    b=10, t=10)
)


fig.show()

Pitch

The pitch of the vector field describes the distance and axis over which a full rotation occurs. It can be defined as the positive eigenvector of $$\Pi_{ij} = \dfrac{1}{4}\epsilon_{ilk}[n_l \partial_k n_j + n_l \partial_j n_k - n_j n_l n_m \partial_m n_k] + \dfrac{1}{4}\epsilon_{jlk}[n_l\partial_k n_i + n_l \partial_i n_k - n_i n_l n_m \partial_m n_k]$$

Still in progress

Vector field

The vector field can be explicitly visualised. It is uniform in the far field with a knotted localised structure near the origin.

x, y, z = symbols('x y z')

# compute the preimages of the z = 0 orientation
vectorx = lambdify((x, y, z), rhopf[0], 'numpy')
vectory = lambdify((x, y, z), rhopf[1], 'numpy')
vectorz = lambdify((x, y, z), rhopf[2], 'numpy')

x, y, z = np.mgrid[-5:5:10j, -5:5:10j, -5:5:10j]

u = vectorx(x, y, z)
v = vectory(x, y, z)
w = vectorz(x, y, z)

fig = plt.figure(figsize=(5,5))
ax = fig.gca(projection='3d')
ax.quiver(x, y, z, u, v, w, length=1)
plt.show()

Free Energy Density Isosurface

We can compare this unusual twist structure to the free energy. In the single constant approximation, the Frank free energy density reduces to $$\dfrac{1}{2}K(\nabla \mathbf{n})^2 + Kq_0\mathbf{n}\cdot(\nabla \times \mathbf{n}),$$ where the constant can now be neglected.

This section is coming soon once I know what I'm doing.

Quaternion Hopf map

This is the version of the Hopf map formulated using quaternions. It connects the map to the rotation group.

The Quaternion_Hopf_Map function returns a vector field by performing the following series of maps:

$\mathbb{R}^3 \xrightarrow{\text{stereo}} S^3 \xrightarrow{qpq^*} \mathbb{H} \simeq S^2$.

In the interest of time, I might abandon this section.

def Quaternion_Hopf_Map(x,y,z):
    pass

Test Functions

These are test functions that should all pass if the code is working properly.

def test_project_Rn_to_Sn():
    
    # R2 to S2
    
    stereo = StereographicProjection(special_point = (0,0,1), sphere_centre=(0,0,0), sphere_radius=1, n=2)
    
    # general point
    vector = (4,5,0)
    output = stereo.project_Rn_to_Sn(vector)
    assert output[0] / (1 - output[-1]) == vector[0]
    
    # origin
    vector = (0,0,0)
    output = stereo.project_Rn_to_Sn(vector)
    assert output == (0,0,-1)
    
    # circle point
    vector = (0,1,0)
    output = stereo.project_Rn_to_Sn(vector)
    assert output == vector
    
    # R3 to S3
    
    stereo = StereographicProjection(special_point = (0,0,0,1), sphere_centre=(0,0,0,0), sphere_radius=1, n=3)

    # general point
    vector = (4,5,1,0)
    output = stereo.project_Rn_to_Sn(vector)
    assert output[0] / (1 - output[-1]) == vector[0]
    
    # origin
    vector = (0,0,0,0)
    output = stereo.project_Rn_to_Sn(vector)
    assert output == (0,0,0,-1)
    
    # circle point
    vector = (0,1,0,0)
    output = stereo.project_Rn_to_Sn(vector)
    assert output == vector
    
    # symbols
    x,y,z = symbols('x y z')
    vector = (x,y,z,0)
    output = stereo.project_Rn_to_Sn(vector)
    assert output[-1] == (x**2 + y**2 + z**2 - 1)/(x**2 + y**2 + z**2 + 1)

test_project_Rn_to_Sn()

def test_project_Sn_to_Rn():
    
    # R2 to S2
    stereo = StereographicProjection(special_point = (0,0,1), sphere_centre=(0,0,0), sphere_radius=1, n=2)
    
    # general point
    vector = (0.5, 0.5, sqrt(2) / 2)
    output = stereo.project_Sn_to_Rn(vector)
    assert np.isclose(output[-1], 0)
    assert np.isclose(output[0] / output[1], float(vector[0] / vector[1]))
    assert np.isclose(1 - float(vector[0]) / output[0], float(vector[-1]))
    
    # origin
    vector = (1,0,0)
    output = stereo.project_Sn_to_Rn(vector)
    assert output == vector
    
    # symbols
    x, y = symbols('x y')
    vector = (x, y, x * y)
    output = stereo.project_Sn_to_Rn(vector)
    assert output[-1] == 0
    assert output[0] / output[1] == vector[0] / vector[1]
    assert 1 - vector[0] / output[0] == vector[-1]
    
    # R3 to S3
    stereo = StereographicProjection(special_point = (0,0,0,1), sphere_centre=(0,0,0,0), sphere_radius=1, n=3)

    # general point
    vector = (0.5, 0.5, 0, sqrt(2) / 2)
    output = stereo.project_Sn_to_Rn(vector)
    assert np.isclose(output[-1], 0)
    assert np.isclose(output[0] / output[1], float(vector[0] / vector[1]))
    assert np.isclose(1 - float(vector[0]) / output[0], float(vector[-1]))
    
    # origin
    vector = (1,0,0,0)
    output = stereo.project_Sn_to_Rn(vector)
    assert output == vector
        
    # symbols
    x, y, z = symbols('x y z')
    vector = (x, y, z, x * y * z)
    output = stereo.project_Sn_to_Rn(vector)
    assert output[-1] == 0
    assert output[0] / output[1] == vector[0] / vector[1]
    assert 1 - vector[0] / output[0] == vector[-1]
    
test_project_Sn_to_Rn()